#########################################################
# Application
# US election
# Replicate Alvarez & Nagler 1995
#########################################################

# load library
library(mlogit)
# read in data
dat = read.table("replication_AN/nes9212r.asc")
varnames <- c("preschc", # choice: 1 Bush 2 Clinton 3 Perot
             "educ", # Respondent's education
             "east", # Region (East)
             "south", # Region (South)
             "west", # Region (West)
             "women", # Gender (Female)
             "respfinp", # Felt personal finance were worse
             "natlec", # Felt national economy was worse
             "bclibdis", # Ideological Distance
             "gblibdis", # Ideological Distance
             "rplibdis", # Ideological Distance
             "dem", # Democrat
             "rep", # Republican
             "respgjob", # Oppose government jobs
             "resphlth", # Oppose government health care
             "respblk", # Oppose government minority assistance
             "respab", # Abortion
             "termlim", # Term limits
             "age1829", # Age: 18-29
             "age3044", # Age: 30-44
             "age4559", # Age: 45-59
             "newvoter", # New or returning voter
             "deficit") # Felt deficit was a major problem

names(dat) <- varnames
dat$choice <- "Bush"
dat$choice[which(dat$preschc==2)] <- "Clinton"
dat$choice[which(dat$preschc==3)] <- "Perot"

names(dat)[9] <- "dist.Clinton"
names(dat)[10] <- "dist.Bush"
names(dat)[11] <- "dist.Perot"

vote_dat <- mlogit.data(dat, varying = c(9:11), shape = "wide", choice = "choice")

vote_dat$bush <- ifelse(vote_dat$alt=="Bush",1,0)
vote_dat$clinton <- ifelse(vote_dat$alt=="Clinton",1,0)
vote_dat$perot <- ifelse(vote_dat$alt=="Perot",1,0)

vote_dat$respfinp.bush <- vote_dat$respfinp * vote_dat$bush
vote_dat$respfinp.clinton <- vote_dat$respfinp * vote_dat$clinton

vote_dat$natlec.bush <- vote_dat$natlec * vote_dat$bush
vote_dat$natlec.clinton <- vote_dat$natlec * vote_dat$clinton

vote_dat$respgjob.bush <- vote_dat$respgjob * vote_dat$bush
vote_dat$respgjob.clinton <- vote_dat$respgjob * vote_dat$clinton

vote_dat$resphlth.bush <- vote_dat$resphlth * vote_dat$bush
vote_dat$resphlth.clinton <- vote_dat$resphlth * vote_dat$clinton

vote_dat$respblk.bush <- vote_dat$respblk * vote_dat$bush
vote_dat$respblk.clinton <- vote_dat$respblk * vote_dat$clinton

vote_dat$respab.bush <- vote_dat$respab * vote_dat$bush
vote_dat$respab.clinton <- vote_dat$respab * vote_dat$clinton

vote_dat$east.bush <- vote_dat$east * vote_dat$bush
vote_dat$east.clinton <- vote_dat$east * vote_dat$clinton

vote_dat$south.bush <- vote_dat$south * vote_dat$bush
vote_dat$south.clinton <- vote_dat$south * vote_dat$clinton

vote_dat$west.bush <- vote_dat$west * vote_dat$bush
vote_dat$west.clinton <- vote_dat$west * vote_dat$clinton

vote_dat$newvoter.bush <- vote_dat$newvoter * vote_dat$bush
vote_dat$newvoter.clinton <- vote_dat$newvoter * vote_dat$clinton

vote_dat$termlim.bush <- vote_dat$termlim * vote_dat$bush
vote_dat$termlim.clinton <- vote_dat$termlim * vote_dat$clinton

vote_dat$deficit.bush <- vote_dat$deficit * vote_dat$bush
vote_dat$deficit.clinton <- vote_dat$deficit * vote_dat$clinton

vote_dat$dem.bush <- vote_dat$dem * vote_dat$bush
vote_dat$dem.clinton <- vote_dat$dem * vote_dat$clinton

vote_dat$rep.bush <- vote_dat$rep * vote_dat$bush
vote_dat$rep.clinton <- vote_dat$rep * vote_dat$clinton

vote_dat$women.bush <- vote_dat$women * vote_dat$bush
vote_dat$women.clinton <- vote_dat$women * vote_dat$clinton

vote_dat$educ.bush <- vote_dat$educ * vote_dat$bush
vote_dat$educ.clinton <- vote_dat$educ * vote_dat$clinton

vote_dat$age1829.bush <- vote_dat$age1829 * vote_dat$bush
vote_dat$age1829.clinton <- vote_dat$age1829 * vote_dat$clinton

vote_dat$age3044.bush <- vote_dat$age3044 * vote_dat$bush
vote_dat$age3044.clinton <- vote_dat$age3044 * vote_dat$clinton

vote_dat$age4559.bush <- vote_dat$age4559 * vote_dat$bush
vote_dat$age4559.clinton <- vote_dat$age4559 * vote_dat$clinton

model <- mlogit(choice ~ dist + respfinp.bush + natlec.bush + respgjob.bush + resphlth.bush + 
                 respblk.bush + respab.bush + east.bush + south.bush + west.bush + newvoter.bush + termlim.bush +
                 deficit.bush + dem.bush + rep.bush + women.bush + educ.bush + age1829.bush + age3044.bush + age4559.bush
               + respfinp.clinton + natlec.clinton + respgjob.clinton + resphlth.clinton + 
                 respblk.clinton + respab.clinton + east.clinton + south.clinton + west.clinton + newvoter.clinton + termlim.clinton +
                 deficit.clinton + dem.clinton + rep.clinton + women.clinton + educ.clinton + age1829.clinton + age3044.clinton + age4559.clinton,
               data = vote_dat,
               reflevel = "Perot",
               method = "bfgs",
               probit=TRUE)

tmp1 <- sqrt(dat$dist.Bush) + 5.32
tmp1[which(tmp1>7)] <- 5.32 - sqrt(dat$dist.Bush[which(tmp1>7)])
tmp2 <- 5.32 - sqrt(dat$dist.Bush)
tmp <- cbind(tmp1,tmp2)
tmp <- cbind(tmp,rep(1,909))
for (i in 1:909){
  if (tmp[i,1] - floor(tmp[i,1]) > 0.1){
    tmp[i,3] <- 2
  }
}
resplib <- rep(0.1,909)
for (i in 1:909){
  resplib[i] <- tmp[i,tmp[i,3] ]
}
resplib <- floor(resplib)

# Correct position:
# Bush: 5.32
# Clinton: 2.98
# Perot: 4.49

# reconstruct dataset
dat$resplib <- resplib
vote_dat2 <- mlogit.data(dat, varying = c(9:11), shape = "wide", choice = "choice")

vote_dat2$bush <- ifelse(vote_dat2$alt=="Bush",1,0)
vote_dat2$clinton <- ifelse(vote_dat2$alt=="Clinton",1,0)
vote_dat2$perot <- ifelse(vote_dat2$alt=="Perot",1,0)

vote_dat2$respfinp.bush <- vote_dat2$respfinp * vote_dat2$bush
vote_dat2$respfinp.clinton <- vote_dat2$respfinp * vote_dat2$clinton

vote_dat2$natlec.bush <- vote_dat2$natlec * vote_dat2$bush
vote_dat2$natlec.clinton <- vote_dat2$natlec * vote_dat2$clinton

vote_dat2$respgjob.bush <- vote_dat2$respgjob * vote_dat2$bush
vote_dat2$respgjob.clinton <- vote_dat2$respgjob * vote_dat2$clinton

vote_dat2$resphlth.bush <- vote_dat2$resphlth * vote_dat2$bush
vote_dat2$resphlth.clinton <- vote_dat2$resphlth * vote_dat2$clinton

vote_dat2$respblk.bush <- vote_dat2$respblk * vote_dat2$bush
vote_dat2$respblk.clinton <- vote_dat2$respblk * vote_dat2$clinton

vote_dat2$respab.bush <- vote_dat2$respab * vote_dat2$bush
vote_dat2$respab.clinton <- vote_dat2$respab * vote_dat2$clinton

vote_dat2$east.bush <- vote_dat2$east * vote_dat2$bush
vote_dat2$east.clinton <- vote_dat2$east * vote_dat2$clinton

vote_dat2$south.bush <- vote_dat2$south * vote_dat2$bush
vote_dat2$south.clinton <- vote_dat2$south * vote_dat2$clinton

vote_dat2$west.bush <- vote_dat2$west * vote_dat2$bush
vote_dat2$west.clinton <- vote_dat2$west * vote_dat2$clinton

vote_dat2$newvoter.bush <- vote_dat2$newvoter * vote_dat2$bush
vote_dat2$newvoter.clinton <- vote_dat2$newvoter * vote_dat2$clinton

vote_dat2$termlim.bush <- vote_dat2$termlim * vote_dat2$bush
vote_dat2$termlim.clinton <- vote_dat2$termlim * vote_dat2$clinton

vote_dat2$deficit.bush <- vote_dat2$deficit * vote_dat2$bush
vote_dat2$deficit.clinton <- vote_dat2$deficit * vote_dat2$clinton

vote_dat2$dem.bush <- vote_dat2$dem * vote_dat2$bush
vote_dat2$dem.clinton <- vote_dat2$dem * vote_dat2$clinton

vote_dat2$rep.bush <- vote_dat2$rep * vote_dat2$bush
vote_dat2$rep.clinton <- vote_dat2$rep * vote_dat2$clinton

vote_dat2$women.bush <- vote_dat2$women * vote_dat2$bush
vote_dat2$women.clinton <- vote_dat2$women * vote_dat2$clinton

vote_dat2$educ.bush <- vote_dat2$educ * vote_dat2$bush
vote_dat2$educ.clinton <- vote_dat2$educ * vote_dat2$clinton

vote_dat2$age1829.bush <- vote_dat2$age1829 * vote_dat2$bush
vote_dat2$age1829.clinton <- vote_dat2$age1829 * vote_dat2$clinton

vote_dat2$age3044.bush <- vote_dat2$age3044 * vote_dat2$bush
vote_dat2$age3044.clinton <- vote_dat2$age3044 * vote_dat2$clinton

vote_dat2$age4559.bush <- vote_dat2$age4559 * vote_dat2$bush
vote_dat2$age4559.clinton <- vote_dat2$age4559 * vote_dat2$clinton

save(vote_dat,file="replication_AN/vote_dat_origin.RData")
save(vote_dat2,file="replication_AN/vote_dat_sim.RData")

############
# simulation
############

pred0 <- predict(model,newdata=vote_dat)
pred0 <- cbind(pred0,rep(0,909))
for (i in 1:909){
  pred0[i,4] <- which.max(pred0[i,1:3])  
}
length(which(pred0[,4]==dat$preschc))
pred00 <- pred0[,4]
pred00[which(pred00==1)] <- "Perot"
pred00[which(pred00==2)] <- "Bush"
pred00[which(pred00==3)] <- "Clinton"

pred000 <- rep(0,909)
pred000[which(pred00=="Perot")] <- 3
pred000[which(pred00=="Bush")] <- 1
pred000[which(pred00=="Clinton")] <- 2
# prediction accuracy:
length(which(pred000==dat$preschc))/909 # 74.1%

preds <- t(as.matrix(table(pred0[,4])/909))
pos_seq <- seq(from = 0.1, to = 7, by = 0.02)

for (i in 1:length(pos_seq)){
  vote_dat2$dist[which(vote_dat2$bush==1)] <- (vote_dat2$resplib[which(vote_dat2$bush==1)] - pos_seq[i])^2
  pred0 <- predict(model,newdata=vote_dat2)
  pred0 <- cbind(pred0,rep(0,909))
  for (i in 1:909){
    pred0[i,4] <- which.max(pred0[i,1:3])  
  }
  preds <- rbind(preds,t(as.matrix(table(pred0[,4])/909)))
}
save(preds,file="replication_AN/pred_voteshare.RData")


